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Resume. — En sous-resultat de l'algorithme de Schoof-Elkics-Atkin pour compter 
le nombre de points d'une courbe elliptique definie sur un corps fini de caracteristique 
p, il existe un algorithme qui, pour £ un nombre premier d'Elkies, calcule des points de 
^-torsion dans une extension de degre I — 1 a l'aide de 0(£ max(f, logg) 2 ) operations 
elementaires a condition que £ ^ p/2. 

Nous combinons ici un algorithme rapide du a Bostan, Morain, Salvy et Schost 
avec l'approche p-adique suivie par Joux et Lercier pour obtenir pour la premiere fois 
un algorithme valide sans limitation sur £ et p et de complexity similaire. 

Abstract. — As a subproduct of the Schoof-Elkies-Atkin algorithm to count points 
on elliptic curves defined over finite fields of characteristic p, there exists an algorithm 
that computes, for £ an Elkies prime, ^-torsion points in an extension of degree £ — 1 
at cost 0(£ maux(£, logq) 2 ) bit operations in the favorable case where I pjl. 

We combine in this work a fast algorithm for computing isogenies due to Bostan, 
Morain, Salvy and Schost with the p-adic approach followed by Joux and Lercier to 
get for the first time an algorithm valid without any limitation on £ and p but of 
similar complexity. 



1. Introduction 

Let K be a finite field with q elements and E be an elliptic curve over K given by 
a plane equation of the form 

(1.1) y 2 + a\xy + a^y — x 3 + a-ix" 1 + a±x + 

where the coefficients a\, a.2, 0,3, 04 and ae are elements of K. For any field L such 
that K C L, we denote by E(L) the set of L-points of E, i.e. the set of solutions 
in L of Equation flTTT]) , plus the additional point at infinity O with homogeneous 
coordinates (0:1:0). The curve E/K. has a structure of commutative algebraic 
group with neutral element O, derived from the secant and tangent rules. Its order 
is equal to q + 1 — t with t € Z such that \t\ ^ 

We are interested in the determination of ^-torsion points of E, that is the set E[£] 
of points P of E(K) such that £P — O for prime integers £, distinct from p. This 
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group is isomorphic to Z/£Z x Ijtl (cf. [171 p. 89]), its cardinal is thus £ 2 . In fact, the 
multiplication by £ is given by a rational transformation of P 2 (K), of degree £ 2 , of the 
form (x : y : z) > (Xi (x, y, z) : Y^(x, y, z) : Zi(x, y, z)) where Xg, Yi and Zi are three 
homogeneous polynomials of degree £ 2 and ^-torsion points are explicitely given by 
Zg(x 1 y, z) = 0. Excluding O, this equation can be easily transformed into an equality 
of the form ft(x) — where fi is an univariate polynomial of degree (£ 2 — l)/2, called 
the £-th division polynomial. 

The improvements by Atkin and Elkies to Schoof 's algorithm for counting points 
on elliptic curve stem from the fact that when the principal ideal (£) splits in the imag- 
inary quadratic field Q(-\/t 2 — 4q), in half the cases thus, there exists two subgroups 
of degree £ in E[i] defined in a degree £ — 1 extension of K. Such an integer £ is called 
an Elkies prime. In this work, we more precisely focus on algorithmic efficient ways 
to compute degree (£ — l)/2 polynomials the roots of which are abscissas of points 
contained in such subgroups. We call these subgroups, and these degree (£ — l)/2 
polynomials over K, £-th Elkies subgroups, and £-th Elkies polynomials. 

Our main result, where we classically denote by <j>i(xi, ... ,Xk) = 0{<f>2{x\, . . . , Xk)) 
functions <fii and 4>2 such that there exists an integer k with (f>i(x\,..., Xk) = 
0((f) 2 (xi, . . . , Xk) log fc (j>i{xi, . . . , x k )), is as follows. 

Theorem 1. — Let E be an elliptic curve defined over a finite field withq elements 
and £ be an Elkies prime, distinct from the characteristic o/K, then there exists an 
algorithm which computes an £-th Elkies polynomial at cost 0(£ max(£, logq) 2 ) bit 
operations and space. 

This problem is closely related to the problem of computing separable isogenies 
of degree £ between two elliptic curves since an application of Velu's formulas [19] 
with starting point such polynomials yields an isogeny. Especially, counting points 
on elliptic curves first raised interest for such computations. But isogenies now play 
a role in numerous other fields, for instance to protect elliptic curve cryptographic 
devices against physical side attacks [18] , to improve Weil descent to calculate elliptic 
discrete logarithms [10] , to decrease the complexity of computing discrete logarithms 
in some family of finite fields [7] , to exhibit normal basis in finite field extensions [6] , 
etc. 

We first recall in Section [5] the complexity of the algorithms known to solve this 
problem. In Section [3] we focus on the fastest algorithm in finite fields of large 
characteristic published so far, due to Bostan, Morain, Salvy and Schost [2j. We 
then show in Section [5] how we can combine this algorithm with the p-adic approach 
introduced by Joux and Lercier in [11] to get a fast algorithm in any finite field and 
we clarify that we need a p-adic precision of only 0(log £/ logp). A detailed example 
is given in Section [5] 
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2. Related work 

For the sake of simplicity, we restrict ourself to finite fields K of characteristic 
larger than three, and to prime integers £ > 2. In this case, an elliptic curve is 
simply given by a plane equation of the form y 2 = x 3 + a^x + ag. Its discriminant, 
always non zero, is equal to = — 16(4et4 3 + 27aQ 2 ) and its j-invariant is equal to 
j E = -12 3 (4a 4 ) 3 /A £ . 

2.1. Naive approach. — £-th Elkies polynomials are factors of the £-th division 
polynomial fg. Therefore, a naive approach consists in computing fi, which can be 
done at cost 0(£ 2 \ogq) elementary operations thanks to a "Square and Multiply" 
method [17], and then in factorizing it with cost 0(£ 1 ' 815x2 log 2 q) |16j . This algo- 
rithm needs a total of (9(£ 4 log 2 q) bit operations. 

2.2. Schoof-Elkies-Atkin framework. — Let tte be the Frobenius endomor- 
phism of E. Its restriction to E[£), seen as a F^-vector space of dimension two, is 
still an endomorphism. When £ is an Elkies prime, its eigenspaces correspond to £-th 
Elkies subgroups C of E[£] and from each C one can construct an isogeny of degree 
£ between E and the elliptic curve E' — E/C, defined over K. 

The following algorithm takes advantage of these facts. 

Step 1 : Compute the modular polynomial of degree £, $>i(X, Y), equation of the 
modular curve Xq (£) . This is a bivariate symmetric polynomial, of degree £ + 1 
in X and Y, whose coefficients are integers of 0(£) bits (cf. |4j). j-invariants of 
^-isogenous elliptic curves are roots of &i(X, Y). 

Step 2 : Compute roots j\ and j2 of &e(X, js). 

Step 3 : Compute a normalized Weierstrass equation for elliptic curves of j- 
invariants ji and j% , and the sum p\ of the abscissas of points in the kernel of the 
isogeny, using the polynomials d^ e /dX, d^f /dY, d 2 &i/dX 2 , d 2 $i/dXdY 1 
d 2 $ e /dY 2 (cf. [15]). 

Step 4 : Compute from each isogenous curve, a £-th Elkies polynomial thanks to 
the kernel of the corresponding isogeny. 

The complexity of the method comes now. 

Step 1 : The modular polynomial $i(X, Y) has 0(£ 2 ) coefficients, each with 
about 0(£) bits. There exists methods to compute this polynomial at cost 
quasi-linear in its size, i.e. 0(£ 3 ) bit operations (cf. [9]). We need to reduce 
this polynomial modulo p, that is 0(£ 3 ) bit operations too. The result is then 
of size 0(£ 2 \ogp) bits. 

Step 2 : With the help of Horner's method, the evaluation of &e(X, Y) at je costs 
0(£ 2 logg) bit operations. In order to compute roots of the resulting degree £+1 
polynomial, we have first to compute its gcd with X q — X, that is 0(£ log 2 q) 
bit operations (cf. |12j ). We obtain a degree 2 polynomial whose roots can then 
be found with negligible cost. 

Step 3 : The computations of the derivatives of $^ and their evaluations can be 
done at cost 0(£ 2 \ogq) bit operations. 
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Step 4 : Here, we have to distinguish several cases. 

— In finite fields of large characteristic, the best algorithm known so far to 
compute isogenies is due to Bostan et al. [2] and takes time 0(1 \ogq) bit 
operations. 

— In finite fields of small but fixed characteristic, the best algorithm known 
is due to Couveignes [5] and needs 0(£ 2 \ogq) bit operations (but the 
contribution of p in the O complexity constant is exponential in logp). 

— In between, that is finite fields of small but non-fixed characteristic, 
the best algorithm is due to Joux and Lercier and needs 0((1 + 
£/p)£ 2 \ogq) bit operations. 

The best total complexity is thus equal to Oil max(^, log q) 2 ), achieved in finite 
fields of large characteristic. But, in finite fields of small characteristic, the complexity 
can be as large as 0(£ 3 \ogq) bit operations when I ^ p. 

This work yields an algorithm of same complexity as in the large characteristic 
case without any limitation on the characteristic or the degree of the base field K. 



3. The large characteristic case 

In order to get an algorithm with good complexity in finite fields of small charac- 
teristic too, we first reformulate the algorithm of Bostan et al. in such a way that 
its extension in the p-adics is more easily studiable. The general strategy is the same 
except that we take into account some specificities of the involved differential equa- 
tion in the resolution. As a result, we obtain a precise and compact algorithm (cf. 
Algorithm [T|) . 

3.1. Differential equation. — In a field K of characteristic larger than three, an 
isogeny between two elliptic curves, E : y 2 = x 3 +a4 x+ae and E' : y 2 — x 3 +a' 4 x+a' 6 , 
can be given by 

(N(x) (N{x) 



D( X y * \D{x) : 

where N and D are unitary polynomials of degree I and i — 1. When c is equal 
to one, the isogeny is said to be normalized. This is in particular the case in the 
Schoof-Elkies-Atkin framework. 

If we now state that the image of a point of E by / is on E', we get the following 
differential equation 

N(x)\' 2 fN(x)\ 3 , (N(x) 



(3 - 1} {x3+aiX + a6 \ D (x)J \D(x)J •-*\D(x). 

This equation can be solved with a Taylor series expansion of N(x)/ D(x) — x in l/x 
for l/x close to 0. The relations obtained thanks to Equation (13. 1|) enable to compute 
by recurrence each coefficient in turn, if the first coefficients are known. It is then 
possible to recover N and D with the help of Berlekamp-Massey's algorithm, or one 



ELKIES SUBGROUPS OF ELLIPTIC CURVE ^-TORSION POINTS 



5 



of its optimized variant. In [2], one takes advantage of a Newton algorithm so that 
the number of coefficients computed at each iteration doubles. 
More precisely, let S be defined by 



At the infinity, N(x)/D(x) has a series expansion of the form x + 0(1). We thus 
have S(x) = x + 0(x 3 ) and this knowledge is finally enough to completely recover 
N(x)/D(x). 

3.2. Resolution. — We consider more generally equations of the form S' 2 = G ■ 
(H o S). In Equation (|3.1|) , we have for instance H(z) = a' 6 z 6 + a' 4 z 4 + 1 and 
G(x) = l/(a6 x 6 + 04 x 4 + 1). We now look for a solution modulo x**, where p, is any 
integer given in input. The way to solve this equation is first to assume that we know 
the solution modulo x d and then, thanks to a Newton iteration, to obtain a solution 
modulo x 2d . After roughly log^i such iterations, one gets the full solution. 

We now present a compact algorithm for this task. Its complexity can be easily 
established, it is equal to 0(p log q) bit operations. Its correctness is slightly more 
difficult to prove and we delay it to Appendix IA1 



Algorithm 1 Solving equation S' 2 = G ■ (H o S), S(0) = a and S'(0) = (3. 

Input: n G N, (a, (3) G K 2 , H 6 K[z], G i K[[x}} 
Output: S £ K[i], a solution of the differential equation modulo x* 1 
d< — 2, U< — 1/(3, J< — 1, V< — 1 
S< — a + f3x+ [(G'{Q) + H'{a)(3 z )/{A(3)} x 2 
while [d < n — 1) do 

U < — U ■ (2 - S' ■ U) mod x d 

V < — (V + J ■ (H o S) ■ (2 - V ■ J)) / 2 mod x d 

J < — J • (2 - V ■ J) mod x d 

S < — S + V ■ J (G ■ {H o S) - S 12 ) (U-J/2)dx mod x^M^+i^) 

d< — 2d 
end while 
return S 



Proposition 3.1. — Let (a,/3) G K 2 where K is a finite field of characteristic p, 
let G be a formal series defined over K, let H be a polynomial defined over K such 
that H(a) = 1 and f3 2 = G(0) ^ 0. Let fi G {1, . . . ,p\, then Algorithm^ computes a 
Taylor series (modulo x^ ) of the solution S of the differential equation 




Equation (|3.1[) becomes 



(o e x 6 + a 4 x 4 + l)S'{x) 2 = 1 + a' 4 S(x) 4 



+ 



a' 6 S(x) 6 . 



S'(x) 2 = G(x)H(S(x)), S{0) 



a, S'(0) =[3. 
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3.3. Full algorithm. — We first compute G(x) = l/(aex 6 + 04 x 4 + 1) mod- 
ulo x thanks to the classical iterative following formula, G\(x) — 1, G2d(x) = 
G d {x) (2 - G d {x) ■ (a 6 x 6 + a 4 x 4 + 1)) mod x 2d . We then apply Algorithm [1] to G(x) 
and H{z) = a' & z 6 + a' 4 z 4 + 1 with fx = U, a = and /3 = 1. 
The obtained solution S is odd, we define from it 

21-1 

T(x) =^2ux\ where Vi € {0, . . . , 21 - 1}, U = s 2i +i. 

We denote by R(x) the inverse of the square of T(x), modulo x 2e , with the same 
inverse formulas as those used for G. We then have 

N(x) fl\ , v x e N(l/x) 

—-4 =xR[-) , i.e. R(x) = , , V A • 

Applying Berlckamp-Massey algorithm [H 1131 [8] or one of its optimized variant [3| 
1141 to i? yields I? and the searched ^-th Elkies polynomial is equal to the square root 
of D. 



4. Extension to any finite field 

To extend the Schoof-Elkies-Atkin framework in any characteristic, the techniques 
developed in [11] give the general idea: to use the p-adics to authorize divisions by the 
characteristic p of the field. These divisions make it possible to use in any finite field 
algorithms primarily designed in large characteristic. There exists one main obstacle 
with this approach. Calculations in the p-adics imply losses of precision at the time 
of divisions by p. It is thus necessary to anticipate a sufficient precision, which results 
in an increase in the size of the handled objects. 

One could hope to perform this lift in the p-adics only in the last stage of the 
algorithm, i.e. for the calculation of the isogeny. It is actually not possible because 
fast algorithms for computing isogenies need normalized models for the isogenous 
curves. 

It is thus necessary to lift in the p-adics from the very beginning of the algorithm. 
It is exactly what is done in with a p-adic precision linear in £. Instead, we 

consider here the techniques of ;2j , and one shows that the necessary p-adic precision 
can be brought back to only 0(log 2 tj logp). The total complexity of the algorithm 
is then similar to the one of the large characteristic case, that is 0(£ max(f, \ogq) 2 ). 

4.1. Lifting curves and isogenies. — One starts by lifting arbitrarily the curve E 
in the p-adics. Any coefficient S4 and clq such that 04 = 04 mod p and olq — mod p 
is appropriate and one works on the elliptic curve E/Q q with model y 2 = x 3 +a4 x+clq. 

The computation of the j-invariant jg of the curve E, of the solutions j\ and j% 
of the equation </>t(x,jE) — 0, as well as Weierstrass models of the corresponding 
curves E\ and E2 , proceeds exactly as in the SEA framework. The curves E\ and E2 
are f-isogenous with the curve E, and the isogenies can be calculated as in the large 
characteristic case. 
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Projection E\ of the curve E\ on the base field K is ^-isogenous with E, and the 
connecting isogeny is the projection on the base field of the isogeny connecting E 
to E\. It is the same for E^. It is thus enough to project the denominators of the 
isogenics on K to identify the required factors of the €-th division polynomial of E. 

4.2. p-adic computations. — From now on, we are interested in the p-adic pre- 
cision of the lift of the elliptic curve E. This precision must be large enough so that 
at the end of the resolution of the differential equation with Algorithm [T] the result 
5* can be reduced in K. 

To this purpose, we need first some definitions. 

Definition. — For any positive integer r, one defines PDiv (p,r) by the largest power 
of p which divides r, PDiv(p, r) = max jfc e N|p fc divides r} . 

We denote by Loss(p, £) the sum Xa<i<iog 2 (4^-i) LpLoss(p, £, i), where 

LpLoss(p, ^, i) = max { PDiv(p, r) 1 2* + 1 ^ r < min(2 l+1 ,4£ - 1)}. 

The following lemma relates the precision needed to the function Loss. 

Lemma 4-1- — Let fj, be the p-adic precision of the coefficients a 4 and clq, then 
when fi > Loss(p, £) the polynomials U , V , J and S computed in Algorithm [7] have 
p-adic integer coefficients. Furthermore the precision of the result S is at least equal 
to {pi — Loss(p, I)). 

Demonstration. — One proves this theorem by recurrence on j, the number of it- 
erations of the loop "while" in Algorithm [TJ We assume that at rank j, ^ j < 
log 2 (4£ — 1), the polynomials U, V, J and S have p-adic integer coefficients and that 
their precision is at least equal to fx — < j< j LpLoss(p, £, i). 
Initialization. In input of the algorithm, we have a = 0, (3 = 1, H(z) — d' 6 z 6 + 
a' 4 z 4 + l and G(i) = l/(ag a; 6 + S4 x 4 + 1). The elements 04 , ag, a' 4 and Sg are integers 
of precision /x and thus G and H are of precision fi too (no division by p occurs in 
the computation of G). The same is true for U, V, J and S. 

Heredity. Let j < log 2 (4f — 1), we suppose the assumption true at rank j — 1. At 
the j th iteration, polynomials U, V and J are updated via additions, multiplications, 
derivations and compositions of the values of U, V, J and S before the entry in the 
loop. All these operations preserve the precision, polynomials U , V and J have thus p- 
adic integer coefficients with precision at least equal to [i — J2i ^ % < j-i LpLoss(p, £, i). 

For S, except the integral operation, the calculations preserve the precision. Coef- 
ficients of the series after the integral operation are inverses of degrees between 2 J + 1 
and min.(2 J+ , 4£ — 1). The largest power of p by which we carry out a division is 
thus LpLoss(p, £, j). The absolute precision of the coefficients of S is thus higher or 
equal to fi — Yli<i<j LpLoss(p, £, i). Furthermore, since this precision is positive, 
each coefficient of S is a lift of the coefficient of the series deduced from the isogeny 
over K, and these coefficients are p-adic integers. □ 
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To minimize the loss of precision, we may use the additional fact that S is odd. 
We thus have to consider only coefficients of odd degree in the algorithm and the loss 
of precision in the loop of the algorithm becomes 

LpLoss '(jp,£,i) = max { PDiv(p,2r + 1) /2 i ~ 1 s$ r min(2 l - 1,2^-1)} . 

Lemma 14.21 yields a clear asymptotic bound on the loss of precision stated in 
Lemma 14. f I 

Lemma 4-2. — We have Loss(p, £) — O (log 2 lj logp) . 

Demonstration. — For all % < log 2 (4£ — 1), LpLoss(p, £, i) is the largest power 
of p which divides a range of integers, at most equal to 2 l+1 , we have therefore 
LpLoss(p, £, i) ^ log p 2 l+1 , and 

Loss(p,£) s? log p 2 (Ei<i<iog 2 («-i) (* + l)) . 

s? log p 2 log 2 (M - 1) (log 2 (M - 1) + 1) , 
^ (log 2 (4^i) + l) 2 /log 2 p. 

□ 

We finally can state our main result. 

Proposition 4-1- — A p-adic precision of 0(log 2 lj logp) is asymptotically enough 
to compute a i-th Elkies polynomial. The total computation needs 0(1 max(£, logq) 2 ) 
bit operations. 

Demonstration. — Computations performed in the Schoof-Elkies-Atkin framework, 
especially calls to Algorithm [1] are realized in the p-adics with precision at most 
0(log 2 e/\ogp). This precision does not modify the O complexities of the large charac- 
teristic case and we still have in the p-adic case a total complexity equal to 0(£ max(£, 
logg) 2 ) bit operations. □ 



5. Experiments 

We have implemented this algorithm in the computer algebra system MAGMA. 
Thanks to it, we were able to observe that the bound on the precision stated in 
Proposition 14. f I is tight. We can illustrate the method with an example too. 

5.1. p-adic precision. — Figure [1] shows the evolution of the precision whenp and 
£ vary. The "The(oretical)" bound mentioned corresponds to Loss(p, £) calculations. 
The "Obs(erved)" bound is what seems necessary at the time of calculations (checked 
on some examples). 

It turns out that the precision observed in practice is near the theoretical bound. 
For many values of £, a gap between the theoretical bound and the observed bound 
appears, but this difference remains quite small. 
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Figure 1. p-adic precisions for p = 5,7, 11 and I < 257. 



5.2. Example. — Let E : y 2 = x 3 + x + 4 be defined over F 5 and £ = 11. 
We first need to compute an upper bound for the 5-adic precision, 

LpLoss(5, 11, 1) = 0, LpLoss(5,ll,2) = 1, LpLoss(5, 11,3) = 1, 
LpLoss(5, 11, 4) = 2, LpLoss(5, 11, 5) = 1 . 

We find Loss(5, 11) = 5 and the 5-adic precision is thus 6. 

A 5-adic lift of the curve is y 2 = x 3 + x + 4. With the help of the classical 
5-th modular polynomial $n, we find that a 11-isogenous curve is given by y 2 = 
.x 3 - 7329.x- 3934 + C>(5 6 ). 

We can now compute the series G(x) modulo x Ai ~ x . 

G(x) = 4374a; 42 + 4298a; 40 - 2331a; 38 - 44 1 7a; 36 + 39 36a; 34 + 35 05a; 32 

+ 228a; 30 - 1041a; 28 - 616a; 26 + 97a; 24 + 236a; 22 + 95a; 20 - 48a; 18 

- 47a; 16 - 12a; 14 + 15a; 12 + 8a; 10 + a; 8 - 4a; 6 - a; 4 + 1 + C>(5 6 ) mod a; 43 . 

A solution of the differential equation based on G{x) and H (z) — a' 6 z 6 + a' 4 z 4 + 1 
is then given modulo x 44 by 

S(x) = - (2 + 0(5)) a; 43 + (2 + 0(5)) a; 41 - (1 + 0(5)) a; 39 + (8 + ( 5 2 )) a; 37 

- (1 + 0(5)) a; 35 + (0(5 2 )) a; 33 + (0 ( 5 2 )) a; 31 - (10 + 0(5 2 )) a; 29 - (7 + 0(5 2 )) a; 27 
- (1 + 0(5 2 )) x 25 + (192 + 0(5 4 )) a; 23 + (125 + 0(5 4 )) x 21 + (293 + 0(5 4 )) a; 19 
+ (4 + 0(5 4 )) a; 17 - (161 + 0(5 4 )) a; 15 - (611 + 0(5 5 )) a; 13 + (211 + 0(5 5 )) a; 11 
- (1494 + C>(5 5 )) a; 9 + (1058 + C>(5 5 )) a; 7 - (733 + C>(5 5 )) x 5 + (0(5 6 )) a; 3 + (l + 0(5 6 )) x , 
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and modulo 5, we find 
T(x) = 3a; 21 + 2x 20 + 4x 19 + 3x 18 + 4x 17 + 3x 15 + 3x 13 + Ax 12 + 2X 11 

+ 3x 9 + 4x a + Ax 7 + 4i 6 + x 5 + x 4 + 3x 3 + 2x 2 + 1 mod x 
We have R(x) = l/T{x) 2 mod x 21 , that is 
R(x) = 2x 20 + 2x 19 + 3x is + x 16 + 2x 15 + 3x 14 + x 13 + 3x 12 + 2x rL 

+ 2x 10 + 2x 8 + 3x 7 + 4x & + 4x 5 + 4x 3 + x 2 + 1 mod x 
The rebuilding of the rational fraction corresponding to R gives 

_ 3X 11 + x 9 + x 8 + x 7 + x 6 + 3x 5 + 2x 4 + 3x 3 + 2x 2 + 2x + 1 
{X > ~ x 10 + x$ + x» + x 7 + 3x 6 + 3x 5 + 3a; 4 + 2 x 3 + x 2 + 2 x + 1 m ° 1 
One reverses the order of the coefficients of the denominator to obtain 

D(x) = x 10 + 2 x s + x 8 + 2x 7 + 3x a + 3x 5 + 3x 4 + x 3 + x 2 + x + 1. 
The ^-th Elkies polynomial is then 

\/D(x) = x 5 + x 4 + x 2 + 3x + 1. 
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Appendice A 
Proof of Proposition 13.11 

Let d be a non-zero even integer, we assume that we know a solution of the differ- 
ential equation modulo x d+1 . We thus have 

(A.l) S' d =G-{Ho S d ) mod x d , S d (Q) = a , S' d (0) = (3. 

Let S 2 d = Sd + A2d be a solution modulo x 2d+1 , with x d+1 dividing A 2d , therefore 
(S d + A' 2d ) 2 = G ■ (Ho (S d + A 2d )) mod x 2d . This yields a linear differential equation 
in A 2d . 

2S' d -A' 2d -G- (H' o S d ) ■A 2d = G-(Ho S d ) - S' d 2 mod x 2d . 
With initial condition A 2d (0) = 0, a solution of this equation is 

/a^ „ 1 f (G ■ (H o S d ) - S' d 2 ) ■ J 2dj 9rf+1 
(A.2) A 2d = -r- / ^ — — dx mod x 2d+1 , 

•J2d J 4& d 

( f G ■ (H' o S d ) , \ , 2d , , 
where J 2 d = exp ^— J dx J mod x + . 

From Eq. ifATT]) , we know that (G ■ (H o S d ) - S' d 2 ) is divisible by x d . Moreover, S' d 
has a non-zero constant coefficient. A factor x d appears then in the integral and it's 
enough to compute J 2 d modulo x d . The inverse of J2d is multiplied by the integral, it 
will thus be multiplied by x d+1 , and it's enough to evaluate this inverse modulo x d . 
The inverse of S' d is needed in the computations of A 2d and J 2d . In A 2d , this inverse is 
multiplied by x d and we then compute a primitive. In J 2d , we compute only modulo 
x d . In both cases, the inverse of S' d modulo x d is enough. This inverse is provided by 
Eq. ||Sj}: 

sZ= G.(Ho Sd ) m0dxd - 
We plug this expression in the computation of J 2d modulo x d , we find 
G-{H>o Sd ) dx = fS' d -(H'oS d ) dxmodxd 



2S' d J 2(HoS d ) 

log(HoS d ) d 

= mod x . 

2 
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We then find the following nice formulas for J 2( j and 1/J 2 d modulo x d , 

J 2d = 7= „ mod x d , -j— = y/H o S d mod x d . 
V n ° Od Jld 




Figure 2. Computation of {7 2 d, Vbd, Jid and S% 



These formulas allow to efficiently compute Sid from Sd and other known quanti- 
ties. 

— From the inverse of S' d , 2 modulo x d / 2 , denoted by Ud, we use a classical Newton 
iteration to compute U 2 d- Since S'd = Sd/2 mod £ d / 2+1 , we have U2d = Ud mod 
a ,d/2 an( j we conl p U t e the coefhcients of U2d thanks to 

U 2d = Ud-{2-S' d - U d ) mod x d . 

— From yJH o Sd/2 modulo x d / 2 , denoted by Vd, and the inverse of Vd modulo 
x d / 2 , denoted by Jd, we compute V 2 d and J 2 d as follows. Getting V 2 d consists 
in computing a solution of v 2 — (H o Sd)(x) = 0. We use 

V M =±(lfc + ^) mods'. 

Jd and Vd are by definition inverses of each other modulo x d l 2 . We obtain the 
inverse W 2 d of Vd modulo x d thanks to Newton formulas too, 

W 2d - Jd ■ (2 - V d ■ J d ) mod x d . 

If we now plug this value in the V 2 d formula, we finally find 

2 V 2d = V d + Jd-(Ho S d ) -(2-V d - J d ) mod x d . 
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Another use of Newton's inversion formula yields J 2 d, 
Jid = Jd-(2-Jd- V 2 d) mod x d . 

Thanks to all these equations, we can compute (U 2 d, V 2 d, Jid) from [Ud, Vd, Jd, Sd)- 
The quantity S 2 d is then obtained from Eq. (IA.2j) . 



We illustrate the corresponding computations in Fig. O 

It remains to obtain initial values, for d — 2. Let 7 be defined by S 2 (x) = a + 
13 x + 7 a; 2 mod x 3 . The series S 2 is solution of the differential equation modulo x 2 
and thus 1 + Af3^ x = G(x) H(a + (3x) mod x 2 . Once derivated, and evaluated at 
x = 0, we obtain 7, and thus the value of S 2l 
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We deduce as well 



U 2 (x) = — mod x , V 2 (x) — 1 mod x and 



J 2 {x) = 1 mod x. 



